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Abstract. In this paper we study the superfluid density of the two component Fermi 
gas in optical lattices with population imbalance. Three different type of phases, the 
BCS-state (Bardeen, Cooper, and Schrieffcr), the FFLO-state (Fulde, Ferrel, Larkin, 
and Ovchinnikov) , and the Sarma state, are considered. We show that the FFLO 
superfluid density differs from the BCS/Sarma superfluid density in an important way. 
Although there are dynamical instabilities in the FFLO phase, when the interaction is 
strong or densities are high, on the weak coupling limit the FFLO phase is found to 
be stable. 
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1. Introduction 

The field of ultra-cold Fermi gases offers a great tool to study many different problems 
of correlated quantum systems. For example, in recent experiments [U [2] [3] H] [5] [6] 
polarized Fermi gases were considered. These systems make it possible to study physics 
in the presence of mismatched Fermi surfaces, and non-BCS type pairing such as that 
appearing in FFLO-states [7J [8] or Sarma-states [9] [10]. These possibilities have been 
considered extensively in condensed-matter, nuclear, and high-energy physics [II] . 

It is also possible to study many different physical problems with close analogs in 
the field of solid state physics using optical lattices. However, unlike solid state systems, 
ultra-cold gases in optical lattices provide a very clean environment. These systems have 
very few imperfections and if one is interested in imperfections, they can be imposed 
easily on the system. Optical lattices are made with lasers, thus the lattice geometry is 
easy to modify [12] [13], HH [15] by changing the properties of the intersecting laser beams. 
For these reasons, in optical lattices one can study various quantum many-body physics 
problems, such as Mott insulators, phase coherence, and superfluidity. The possibility 
of a superfluid alkali atom Fermi gas in an optical lattice has been recently studied both 
theoretically PSI HZl QSl IISl EDI ED 122] , as well as experimentally [23]. 

The pairing predicted by the BCS theory does not always mean that the gas is 
a superfluid [21]. When a linear phase is imposed on the order parameter, the phase 
gradient corresponds to superfluid velocity, and a part of the gas, which has superfluid 
velocity is called superfluid, and the coefficient of the inertia of this moving part is 
called superfluid density [25], [26] [27] . However, if the gas is normal the order parameter 
vanishes and thus the phase shift does not affect to the normal gas. If the superfluid 
density is positive in all directions the gas is a superfluid. Negative superfluid density 
implies dynamical instability of the gas [28] [29] . 

There are a fair number of studies about superfluid density of Fermi gas in the free 
space as well as in a trap [30] [31] [32] [33] , and few studies about superfluid density in 
optical lattices [34] [29]. However, these papers have focused on population balanced 
cases i.e. the densities of the components are same. In this paper we study the 
superfluid density of an ultra-cold two component Fermi gas in optical lattices at finite 
temperatures and with finite polarization. We are motivated by the fact by calculating 
the superfluid density we can draw more conclusions on whether the gas is actually 
superfluid or not. Furthermore, we investigate if the superfluid density can be used 
as an indicator of different superfluid phases. The phases we consider are the BCS 
phase, i.e. the densities of the component are equal, the one mode FFLO phase, i.e. 
the densities are not equal and the phase of the pairing gap modulates as a function 
of position, and Sarma phase i.e. the densities are not equal and the phase of the 
pairing gap is constant. We show that there are qualitatively difference between the 
BCS/Sarma superfluid density and the one mode FFLO superfluid density. We also 
study the stability of different phases. We also demonstrate that there can be dynamical 
instabilities in the one mode FFLO phase. 
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This paper is organized as follows. In Sec. [2] we discuss the physical system and 
present the Hamiltonian of the system. In Subsec. 12. li the mean-field approximation and 
the ansatzes we consider are presented. In Sec. [3] we determine the superfluid density 
tensor. In Sec. H] we present the numerical results and we end with some concluding 
remarks in Sec. [5] 



2. Lowest band fermions in an optical lattice 

We consider a two component Fermi gas, whose components are two different hyperfine 
states of the same isotope, and we call them f-state and j-state. The Hamiltonian of 
the system is given by 

H=J2[ rfr ^(r)(-^ + K(r)-/, CT )i ff (r) (1) 

+ g J J drdr'^{(r)^l(r')5(r-r')^ 2 (r')^ 1 (r) 

where h~ = h/2n, h is Planck's constant, ^(r) and ^(r) are the fermionic creation- 
and annihilation field operators of the component a, [i a is the chemical potential of the 
component a, and 5(r — r') is Dirac's 5-function. The interaction strength is related to 
the s-wave scattering length a s through 

4nh 2 a s 

9 = • 

m 

The lattice potential has a cubic structure and is given by 
V a (r) = E r ^ s *,a sm 2 (kx a ), 

a 

where E r = h 2 k 2 /2m is the recoil energy (k = (vr/c?) where d is a lattice constant), and 
s ata is the lattice depth in the a-direction for the component a. 

When we assume that only the lowest band is occupied the field operators can be 
expanded by using the localized Wannier functions in the following way 

i(r) = ^w ff ,,(r)c ff)i (2) 

% 

i 

where w CTj j(r) is a Wannier function, which is localized at a lattice point i, and c a>i is an 
annihilation operator. 

The lowest band Hubbard model is valid, when the lattice is deep enough. In other 
words, it is valid, when the Wannier functions decay within a single lattice constant, 
and when the effective interaction between atoms is much smaller than the bandgap 
between bands. These two conditions imply that one has to take into account only the 
on-site interactions. In this case overlap integrals of the kinetic energy operator between 
the next nearest neighbor are small compared to overlap integrals of the kinetic energy 
operator between the nearest neighbor [35], and consequently, one needs to take into 
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account only hopping between the nearest neighbours. Then the lowest band Hubbard 
Hamiltonian is given by 

h = - J2 [ J ^Y1 +J ^y J2 £ -A; ( 4 ) 

- X)(^t4 A< + A*) + U J2 4Aaa 



',11 



where (i,j) a means sum over the nearest neighbours in the a-direction, and hopping 
strength is defined by 

f ( h 2 V 2 \ 

The coupling strength of the above Hubbard model is given by 
U = gJ dr\w u (v)\ 2 \ Wll (r)\ 2 - 

2.1. Mean- field approximation 

Because the interaction term U ^ c! ^cj jC^c-^ is hard to handle, we can approximate 
it by using mean-field approximation. Under this approximation the interaction part 
becomes 

Y^^YAA, + A *(*)qAi - 

i 

where the pairing gap A(i) = [/(cj^c-^). With respect to position dependence, we only 
consider the case, in which only the phase of the gap A(z) = |A|e iqRi (Rj is a lattice 
vector) can depend on the lattice site i. This ansatz is called the one mode FFLO (Fulde, 
Ferrell, Larkin, and Ovchinnikov) [3, [8]. This ansatz includes also the BCS-ansatz and 
the Sarma-state ansatz as special cases. In these phases the momentum q is simply 
zero. Under this mean. field approximation the Hamiltonian is given by 

C=1U \ (i,j)x (i,j)y (i,j)zj 

- I>t4a>* + »A/iJ + E i A i e<q ' Ri 4,<4,< + i A i e ~ iq ' Ria i a- 

i i 

By minimizing the free energy Fq = Qq + /ifiVf + /xj,./Vj_ (Qq is the grand canonical 
potential of the mean- field Hamiltonian) with respect |A| and q, and solving the 
number equations dF /dfi a = simultaneously, one finds the pairing gap, the chemical 
potentials, and the momentum q as functions of the temperature and the particle 
numbers N a . 
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3. Superfluid density 

Landau's two component model for a superfluid gas says that the superfluid gas consists 
two components the normal component and the superfluid component. In the free space 
the superfluid density is density of the superfluid component. The case is a little bit 
different in lattices thus the energy difference between the twisted system and the system 
without the twisting phase can depend on the direction of the phase shift. This is related 
the fact that the effective masses can depend on the direction (the effective masses are 
related to the hopping strengths). Thus in the lattice superfluid density is more like a 
tensor than a scalar. 

In the free space when the superfluid gas flows the kinetic energy of the gas is given 

by 

Ek = \ I rfr ^( r K( r ) 2 ' 

where p s (r) is the superfluid density and v s (r) is the superfluid velocity. The superfluid 
velocity is defined by 

v,(r) = A V «Hr), 

where <p(r) is the phase of the order parameter. 

To impose the linear phase variation to the order parameter, one can impose a 
linear phase variation © Rj = (Q x /(M x d), Q y /(M y d), Q z /(M z d)) ■ R; [261 E3 EE] to the 
Hamiltonian. Here M a indicates the number of lattice sites in direction alpha, i.e., the 
total volume the lattice V = M x M y M z d 3 This variation corresponds small (Q a is small) 
superfluid velocity, which is given by 

h_ 

2m ' 

Thus this imposed phase gradient gives the system a kinetic energy, which corresponds 
to the free energy difference F® — F , where F® is the free energy within the the phase 
variation and Fq is the free energy without the phase variation. The superfluid fraction 
in this case can be determined as [26J, [27] 

(6) 

0=0 



—(e x /(M x d), e y /(M y d), e z /{M z d)). 



.. IF & -F 1 d 2 F & 
p aa > = lim 



e^o N J x Q a Q Q/ NJ X dQ a dQ c 
where N is the total number of particles, and J x = (J-[ )X + J{,x)/^- This J x corresponds 
the effective mass 

h 2 

me " = 2jjr 

When one takes the limit — ► the temperature, and the number of particles (and of 
course the interaction strength) should be kept as constants. 
Of course, one could define the superfluid fraction as 
1 F & - F 
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where J a ^ = (Jj- )Q + J^ a )/2, but in this case it is a little bit unclear how to define the 
off-diagonal elements. This definition is not directly connected to the energy difference 
between the twisted system and the untwisted system. 

By imposing the phase variation to the order parameter the one finds 

h® = - E J ^ E +J ^y E E ( 8 ) 

CT =T4 \ (i,j)x {i,j) y 

- 5>t4 A* + Wit*) + E I A| C <fc +M >*4 A + |A|e~ i(q+20),Ri c ii iC T:i . 

i i 

One can make the unitary transformation [26j (27J 36J, which corresponds the local 
transformation c a ^ — > c^e 1 & and the twisted Hamiltonian becomes 

H® = - E(^T c tA< + Vl c \A^ ( 9 ) 

i 

+ i A i eiq ' Ri 4/l + i A i e " q ' Ria iA* 



-E[E( J ^( c<e " /JM '*A^+ e 

<t n 

+ J (e i0z/Mz d^ c + e~ i0z/Mz c f c 



The twist angles G Q have to be sufficiently small to avoid effects other than the collective 
flow of the superfluid component. Since Q a /M a is small, we can expand up to second 
order 

p ±ie a /M a _ 1 i - U a U « 

~ M a 2M 2 " 
In this way we can write the twisted Hamiltonian as 

r 6 2 

H & = H + H' ^H + J2 E J <^(4,A,n+a + cl n c^ n _ di ) + (10) 

a x n 

Q 2 Q2 

2ikf 2 ^ ' J o-,y($cr,nC<j,n+dy + ^cr,n^cr,n-dy) + 2M 2 ^ ' ^ Cr,Z (^n^n+cfl "I" ^a-,n^<y,n-dz) 
" V n ~" 



X -n » 



In the momentum space this formula becomes 

e 2 



i^^o + X: £ 2 J CT , a cos(M)4,A,k (11) 



k 



+ E 2,/<T ' a sin ( fc a rf )4,kC<7,k 

k 

= # + T + J, 
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where 

f = E fy§2 E 2 ™(M<iA,.k| (12) 



and 



L2M2 k 



^ = E [^T E 2J ^ a sin(M)4,k g <x,k • (13) 



These two terms are proportional to the current, and number operator of the system, 
respectively. These two terms commute with the number operators, thus when one 
imposes the perturbation in the system the number of particles is conserved. Using the 
same method, what is used in reference [30], can be shown that in the mean- field level 



lim *m. = lim (<*m. , (14) 



|©Ho dO a d® a > |e|-»o \dQ a dO a// AiAtm 

Therefore we can use the grand canonical potential. 

The grand potential of the system can be written with a perturbing phase gradient 
as a series [37] 

n & = n - V i— ^. / dn ■ ■ ■ / dr n (t t h'( Ti ) ■ ■ ■ H'(r n )) , (15) 

n=1 P ntl Jo Jo 

where (3 = 1/ksT, k^ is Boltzmann's constant, Q is the grand canonical potential in 
the absence of the perturbation. The symbol T T orders the operators such a way that r 
decreases from left to right. The brackets (. . .)o = Tr[exp(— (3(Hq — — /xj_iVj)) . . .] 
mean the thermodynamic average evaluated in the equilibrium state of the unperturbed 
system at temperature T. Because all the twisted angles a are small, we can safely 
ignore terms of higher order than 0^,. With this approximation the grand potential 
becomes 

^^^0+—^ dT ( T r f ( T ))o-^ 2 J Q J Q drdr'(T T J(r)J(r')) ,(16) 

The first perturbation term is easy to calculate and it is given by 

1 r 2 i 

— / dr (T T f(r)) = E[^#E 2 J ^ ™<k a d)N aM \ , (17) 

a, a a k 

where N a ^ is number of particles of the cr-component in the momentum state k. Using 
Wick's theorem the second perturbation term is given by 

drdr' (T t J(t)J(t')) 



2(3h 2 Jo Jo 

I rph r/3h 0/ 

l I dTdT '[ E ]^E 4J -^'-' sin (^) sin (^ rf ) 

cr,<j' ,a,a' k,k' 

((1 - 5 aa ,)5 k _ k / +q F(k, r, T')F j (k', r, r) - <W<WG>(k, r, r')Gv(k', r, t 



x 
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G (T (k,r,r') = -(T r c CT , k (r)ci ik (r'))o = ^= £ e'^^G^ u n ) (19) 



rt=— oo 

oo 



F(k,r,r') = -(TAk +q (r)c i ,_ k+q (r'))o = -^ £ V(k,^ n ). 

v n=— oo 

The fermionic Matsubara frequencies are given by 
_ n{2n + 1) 

where n is an integer. The individual Green's functions of the mean-field theory can be 
written as 

G T (k, u, n ) = - K 'f + - '^'f (20) 
G^u n )= K 'f + '^f (21) 

2/?U; n — -&-,k,q «"^n + -^+,k, q 

£ {k.,uj n ) - — — — , {22) 

lfiw n + k,q ifUVn - -&+,k,q 

where the quasiparticle dispersions and amplitudes are given by 



£ ± kq = ^1 |^T,k + q/2 + 6W q /^ 2 + | A |2 ± ^T,k + q/2 ~ ei.lc+q/2 ^ 



Pk. 



= I + 6 T,k+q/2 + q-k+q/2 \ ^ 
2 V ^+,k,q + k,q / 

K,q| 2 = 1 - |Mk,q| 2 , (25) 

and furthermore e CTi k = z~2 a 2J<r, a (l — cos(k a d)) — /v By calculating the integrals and 
the Matsubara summations one finds a lengthy formula 

a, a' a Q k 

/ 2 , ,2, |2 /(^ + .k,q)+/(^-,k, q )-l 
X U|Mk,q| |Vk, q | — 

V -^+,k,q + -^-,k, q 

/3| Mk ,q| 4 /(^+,k,q)(l " /(£ + ,k,q)) - /5|^k,q| 4 /(^-,k,q)(l " /(^-,k,q))) 

+ 2J La J La , sin((k a - q a /2)d) sm((k a > - q a > /2)d) 

v f 9 U l 2 U, |2 /(^+,k,q)+/(^-,k,q)-l 

x (2|u k>q | |uk, q | — — 

V -fr+,k,q + -^-.k.q 

/3| Mk ,q| 4 /(^-,k,q)(l " /(£-,k,q)) - P\v^f{E + ^) (1 - /(£ + , k , q ))) 

+ 4J TjQ J ijQ / sin((/c a + q a /2)d) sin((A; Q / - q a >/2)d) 



2|„. |2 



X «k,q ^k 



/2(/(£ + , k , q )-/(£_, k , q )) 



qi rMi l p 771 

^ ^+,k,q - ^-,k,q 
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2-E , +, k , q 2EL jkiq 

where f(E) is the Fermi-Dirac distribution. In order to handle the limit lim k /^ k , we 
have implicitly assumed in equation fl26l) that the lattice is large, i.e, M a ^> 1. 
Now the twisted grand canonical potential can be written as 

*"* + S>-TO < 27) 

a,a' 

In all the cases we consider SQ aa r = when a ^ a'. The off-diagonal terms can be non- 
zero only when the single particle dispersion couples the different directions together. 
In our case where the directions are independent, the k sums in different directions 
are independent, and because sine is an antisymmetric function these sums vanishes. 
Therefore the grand potential becomes 

J]@«0 Q + ^5Q aa ^. (28) 

a a 

We can now determine the components of the dimensionless superfiuid fraction tensor 

as 

paa ' = j x (N, + Nl y 

As a formula these components are given by 

Paa = a? ^ a ' a C0S {k a d)N ajk (29) 

k,cr 



+ -^J2[ J l^ 2 ((k a + q a /2)d) 

x (2\u, l 2 l„ ^(E^+fjE^)-! 

x I Z Fk,q| |^k,q| rp , p 



+ 2Jf Q sin 2 ((fc a - g Q /2)d) (2|u kjq |> ki , 



,k,q T -^-,k,q 

/3|«k,q| 4 /(^ + ,k,q)(l " /(£+,k, q )) " /5|^k,q| 4 /(^-,k,q)(l " /(^,k,q)) 

, /(E +|k , q )+/(E_, k , q )-l 
-^+,k,q + -E-,k,q 

/5|«k,q| 4 /(^,k,q)(l - /(£-,k, q )) - /9|^k,q| 4 /(^ + ,k,q)(l - /(^ + ,k,q)) 

+ AJ^ a J i:a i sin((fc a + q a /2)d) sin((fc a - q a /2)d) 

, |2 , |2 / 2(/(g +|k , q )-/(£;-, k , q )) 

X |M k) q| |Vk,q| I p p 

V ^+,k,q - -^-,k,q 

2/(^+,k,q) - 1 2/(S_ 4e , q ) - 1 



2-E+, k , q 2S_ kq 

Paa' = 0, 

where N = + N± and J CT)Q! = J a , a /J x . When the momentum q is non-zero the 
superfiuid density describes the one mode FFLO phase, and when the momentum q = 
the superfiuid density describes the BCS or the Sarma phase. p aa is a superfiuid fraction, 
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it describes what fraction of the atoms is in the superfluid state. In the limit where 
J], a — Ji,a' — J for all a, a' and — N^, the chemical potentials are same and q = 0. 
In this limit superfluid density tensor can be simplified as 



k,cr k 
Paa' = (31) 



where 



^k=W( eT ' k+ 2 £i '- k V + |A| 



In the continuum limit i.e. d —>■ 0, where Jd 2 remains constant, the superfluid 
fraction becomes 



Paa = P = 1 + 77 > k 

m. n v ^— ' 



m a nV ^ z <9£ k ' 

k 

where the effective mass m a = h 2 /(2Jd 2 ), n = n-f + is the total number density and 
V is the size of the system. This continuum result is precisely the Landau's formula for 
the BCS superfluid fraction [38] . 



4. Results 



4-1- BCS/Sarma-phase results 

In the continuum the formula of the BCS superfluid fraction is the well known Landau's 
formula. Figured] (a) shows the BCS superfluid fraction as a function of the temperature, 
in the uniform case (the external potential is zero), without the lattice. In figured] (b) 
we show the BCS superfluid fraction divided by | A | 2 as a a function of the temperature. 
One sees from figure [1] (a) that the BCS superfluid fraction is one at zero temperature 
and from figure [H (b) that the superfluid fraction is almost proportional to |A| 2 , 
but the temperature dependence of the superfluid fraction differs somewhat from the 
temperature dependence of | A| 2 . The standard BCS result is that the superfluid fraction 
is proportional to |A| 2 in limit T — > T c [39]. When we calculated the gap, we used the 
renormalization, i.e., we have cancelled out the divergent part of the gap equation, in 
reference [39] is used the cutoff energy. This is the cause of little difference between the 
results. 

When all the hopping strengths are the same i.e. J = J^ a = J^ a ', and filling 
fractions n a = N a /M x M y M z are same and the average filling fraction n av = (nj + n^)/2 
goes to zero, the superfluid fraction should go to one at zero temperature. For low filling 
fractions the the atoms occupy only the lowest momentum states and the single particle 
dispersions e CTj k = 27^^(1 — cos(k a d)) — fi a approaches to form J{kd) 2 — fi a , which 
is the free space dispersion. Therefore the results on the limit of low filling fractions 
become the free space result and we know from the Landau's formula that the BCS 
superfluid fraction in the free space is one at zero temperature (see figured] (a)) . 
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In figure [2] (a) we show the BCS superfluid fraction as a function of the total filling 
fraction in the lattice. As we can see, the BCS superfluid fraction approaches one, when 
the filling fractions become small. In the free space the superfluid fraction is one at zero 
temperature, but in the lattice this no longer holds. In figure |2] (b) we show p{n^ + 
as a function of the total filling fraction and as one can see from it, the BCS superfluid 
density, which is not divided by the total filling fraction is symmetric with respect to 
half filling. The cause of this is the particle-hole symmetry of the lattice model. The 
particle-hole symmetry of the lattice model can be seen also from figures [2] (c) and [2] 
(d). In figures [2] (c) and[2] (d) we show A/ J and | A| 2 / J 2 , respectively, as functions of of 
the total filling fraction. By comparing figures [2] (b)-(d) one can notice that p{n^ + nj_) 
is not directly proportional to A or |A| 2 . 

Figure [3] (a) shows the BCS superfluid fraction as a function of the temperature, 
with three different average filling fractions. When one compares this figure to figure [1] 
(a), one notices that the lattice result is very similar compared to the free space result. 
The superfluid fraction goes to zero, when the paring gap goes to zero as expected. The 
figure is consistent with figure [2] (a), in that the total filling fraction increases when the 
superfluid fraction decreases. In figure [3] (b) we show the BCS superfluid fraction as a 
function of the dimensionless coupling strength —U / J at zero temperature. As we can 
see, when —U/J increases the superfluid fraction decreases contrary to the free space 
result where superfluid fraction is a constant as a function of the coupling strength at 
zero temperature. Of course at finite temperatures the free space superfluid fraction 
decreases with increasing coupling strength. This is due to the fact as —U/J increases 
the movement of the atoms decreases in the lattice. In other words when —U/J increases 
then Cooper pairs become smaller and the atoms are better localized. Of course in the 
limit —U/J — > oo the mean-field theory fails and in large values of —U/J the results 
are not reliable. Figure [3] (c) shows the BCS superfluid fraction divided by |A| 2 as a 
function of the temperature, with three different average filling fractions and we see that 
the BCS superfluid fraction is almost proportional to |A| 2 (compare to figure [1] (b)) i.e 
p = c(T, n-f = ri|)|A| 2 , where c(T,n^ = njj depends weakly on the temperature. We 
show the pairing gap as a function of —U /J at zero temperature in figure [3] (d). As one 
can notice, the BCS superfluid fraction does not behave anything like the pairing gap 
as a function of —U/ J at zero temperature. 

Since it is experimentally easy to study anisotropic lattices and use these to explore 
dimensional crossovers, let us explore superfluid fractions in anisotropic lattices. We 
show in figure H] (a) the BCS superfluid fraction as a function of the temperature in 
the case, where the hopping strengths are different in the different directions. As one 
can see from the figure the components of the superfluid fraction tensor are different 
in different directions. Furthermore, we see that the superfluid fraction is larger in the 
direction of large hopping strength. This is easy to understand, and is due to the free 
energy difference in a-direction being proportional to J a as we see from equations fl28l) 
and ( 1291) . If the superfluid fraction was defined as in equation ([7]), the components of 
the superfluid fraction would be equal in every direction. In other words the effective 
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masses are different in different directions, figure H] (c) shows that in the in the case, 
where the hopping strengths are different in the different directions, p aa /\A\ 2 is almost 
a constant function of the temperature. 

Figure H] (b) shows the BCS superfluid fraction as a function of the temperature in 
the case, where the hopping strengths are different for the components (but J^+Jj, is kept 
constant). As it is clear from figure H] (b), when the ratio J\/ J\ increases the superfluid 
fraction decreases, but the pairing gap, however, does not decrease as a function of 
J\/J\ at zero temperature (it does as a function of the temperature), but it remains 
almost a constant (at zero temperature it is a constant). This can be seen from figure H] 
(d), in which is shown the pairing gaps as as functions of the temperature in the cases, 
where the hopping strengths are different for the components. One can also notice that 
when Jf/Jj. — > oo then p — > ( in this case the critical temperature also goes to zero). 
When J]/Ji increases but 2U/(J^ + Jj_) remains constant, the lattice becomes deeper 
and deeper for j-component. Thus the atoms of the component j are more localized 
and the atoms do not move easily. 

Figure [5] (a) shows the superfluid fraction as a function of polarization P = 
(nj— njj / '(n-r+nj,), at a constant temperature (ksT/J = 0.75). When the polarization is 
larger than zero the state is called the Sarma state. One sees, that when the polarization 
increases the superfluid fraction decreases. Figure [5] (b) shows the superfluid fraction 
divided by |A| 2 as a function of the polarization at a constant temperature. When 
P is about 0.35 the superfluid fraction divided by |A| 2 drops to zero suddenly. This 
happens because the pairing gap vanishes at this polarization. From figure [5] (b) is seen 
that the superfluid fraction at a constant temperature is almost proportional to |A| 2 . 
Figured] (c) shows the Sarma state superfluid fraction as a function of the temperature, 
at a constant polarization (P = 0.10). Figure [5] (d) shows the Sarma state superfluid 
fraction divided by |A| 2 (p/|A| 2 ). We notice from figures [5] (b) and (d) that p/|A| 2 is 
almost a constant as a function of the temperature and polarization i.e. p = c(T, P) | A| 2 , 
where c(T, P) depends weakly on the temperature and the polarization. 

4-2. Results for FFLO-phase 

When the lattice is cubic and all the hopping strengths are equal the one mode FFLO 
superfluid fraction is symmetric like the BCS/Sarma superfluid fraction. This due the 
following fact, when one varies the q, it turns out that the free-energy is minimized, 
when the q lies along side the axis i.e. q = q a x ai but the system does not favour any of 
these axis. 

The BCS - and the Sarma state superfluid fractions are almost directly proportional 
to | A 2 1 i.e p ~ |A 2 |. This changes for the one mode FFLO superfluid fraction as we will 
now demonstrate. We show in figure [6] (a) p xx -component of the FFLO superfluid 
fraction as a function of the total filling fraction, with three different interactions 
U/J = —4.0,-5.0,-6.0, at zero temperature, figure [6] (b) shows p xx (n^ + n\) as a 
function of the total filling fraction, when U/J = —6.0, at zero temperature. The data 
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which is used in figure E] (a) and (b) have been calculated just above the Clogston 
limit [40] . The Clogston limit is the limit for the chemical potential difference below 
which one can find only the BCS type solution when one minimizes the grand potential. 
The numerical value of this limit is roughly 5 fx ~ \/2Ao, where Ao is the pairing gap at 
zero temperature when 5ft = 0. As one can see from figure [6] (a) when \U/J\ increases 
the superfluid fraction decreases and eventually becomes negative. One can also see, 
that as the total filling fraction increases the superfluid fraction decreases but does 
not necessarily become negative. Negative superfluid fraction implies that the gas is 
unstable. If superfluid fraction is negative the state around which we expanded is not 
the ground state of the system, and thus the state cannot be stable. In the case where 
superfluid fraction is negative, the real minimum of the free energy is beyond our ansatz. 
The real minimum of the system may be a some kind of phase separation or modulating 
gap state (like the two-mode FFLO state) [4T1 |22j |33] . 

The superfluid density i.e. p xx (n^ is symmetric as a function of the total filling 
fraction, with the value + nj_ = 1.0, this is shown in figure [6] (b). This is again due 
to the particle-hole symmetry of the lattice model. 

Figure [7J (a) shows the superfluid fraction as a function of polarization at a constant 
temperature {k B T / J « 0.23), and figure [7] (b) shows the superfluid fraction divided by 
|A| 2 as a function of polarization. There is a second order phase transition between the 
Sarma/BCS phase and the FFLO phase, when polarization P « 0.15, this can be seen 
from the bend in the superfluid fraction, clearly. We see also from figure [7J (b) that the 
FFLO superfluid fraction is not proportional to | A| 2 . When one writes p = c(T, P)\ A| 2 , 
the polarization dependence of c(T, P) is very different in the FFLO phase compared to 
that in the Sarma phase. 




0.01 0.02 0.03 0.04 0.05 0.06 0.01 0.02 0.03 0.04 0.05 0.06 



T7T F ~ " T7T F 

Figure 1. Figure (a) we show the BCS superfluid fraction as a function of the 
temperature in a free space (absence of the lattice). Figure (b) shows the BCS 
superfluid fraction divided by |A| 2 as a function of temperature. The interaction 
strength kpa s = —0.30, where /cf is the Fermi wave vector. 

In figures M and [9] we present phase diagrams with four different interaction 
strengths, when the average filling fractions are n av = 0.2 and n av = 0.5, respectively. 
The figures show that there can be unstable regions in the FFLO phase. In the 
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Figure 2. In figure (a) we show the BCS superfluid fraction p = p xx — p yy — p zz as 
a function of the total filling fraction n-f + n± = 2ri| = 2n^ at zero temperature. In 
figure (b) is shown p(n-f +nj.) as a function of the total filling fraction. Figures (c) and 
(d) show A/ J and |A| 2 /J 2 , respectively as functions of the total filling fraction. All 
the hopping strengths are same i.e. J = J|. Q = J[. a >, and —U/J = 6.0. 



unstable FFLO region the superfluid density is negative, which implies that the phase 
is dynamically unstable. However, these instabilities occur only when the interaction 
is relative high. By comparing the two figures one notices that the unstable regions 
increase when the density increases. Thus when the density decreases the unstable 
regions eventually disappear, and this implies that the FFLO phase region is stable 
without the lattice. Due to the particle-hole symmetry average filling fraction 0.5 is 
the worst case for the stability, but even when the average density is 0.5, the one 
mode FFLO is stable when \U/J\ < 4.4 for any polarization. The induced interactions 
might make the FFLO phase more stable, since these induced interactions make the 
effective interaction weaker [Sj. References [221 HH] considered also the possibility of 
phase separation. However, this work we do not consider phase separation, because 
when there is phase separation between the BCS phase and the normal gas, there is a 
superfluid gas in a part of the lattice and normal gas in another part of the lattice. In 
this case the superfluid density is the standard BCS superfluid density and the gas is 
stable. It should be notice that the region of the phase separation does not include the 
unstable regions of the FFLO phase [22]. 
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Figure 3. Figure (a) shows the BCS superfluid fraction p = p xx = p yy = p zz 
as a function of the temperature, with three different average filling fractions n av = 
+ 7ij)/2. Figure (b) shows the BCS superfluid fraction as a function of a 
dimensionless coupling strength — U/J at zero temperature. In figure (c) we show 
the BCS superfluid fraction divided by |A| 2 as a function of the temperature,with 
three different average filling fractions. Figure (d) shows the pairing gap as a function 
of — U/J at zero temperature. In the all figures all the hopping strengths are same 
i.e. J = Jf iQ = Jx, Q '. In figures (a) and (c) — U/J = 6.0, and in figures (b) and (d) 
= ni = 0.5. 



5. Conclusions 

In this paper we have presented, at the mean-field level, the superfluid density of the 
two component Fermi gas in an optical lattice. We have shown that the BCS superfluid 
density in optical lattices differs from the free space results. We have also shown that 
the one-mode FFLO superfluid density differs crucially from the BCS/Sarma superfluid 
density. In the BCS/Sarma phase the gas is always a stable superfluid, but in the 
FFLO phase dynamical instabilities can appear. However, the FFLO phase is stable, 
when \U/J\ < 4.0 i.e, on the BCS limit. Even in the BCS phase the superfluid density 
can be different in different directions depending on the lattice structure. 

Although the one mode FFLO phase is stable when \U/J\ < 4.0, it is unlikely that 
the one mode FFLO is the free energy minimum the system; calculations made by using 
the Bogoliubov-de Gennes equations show that the spatially modulating gap might be 
more favorable [4"Ij I4"2| H5]. Such states are also, quite likely, more stable than the one 
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Figure 4. [Colour online] In figure (a) we show the BCS superfluid fraction as a 
function of the temperature in the case, where the hopping strengths are different 
in the different directions. In figure (b) we show the BCS superfluid fraction as a 
function of the temperature. Figure (c) shows the BCS superfluid fraction divided 
by |A| 2 as a function of the temperature. In figure (d) we show the pairing gaps 
as functions of the temperature. In figures (a) and (c) J x = («7f, x + Ji,x)/%> 
J hv = Jl,v = °- 8 ^T,^ = 0- 8 ^i,:r, J],z = J[,z = 0.5J T:a; = 0.5J; >a; , n T = 7?.^ = 0.2, 
and —U/J^ x = 6.0. In figures (b) and (d) J a , a = J<r',a' = J<r, —2U/(Jf + Jj_) = 6.0, 
and rij = ni = 0.5 



mode FFLO phase. 

The methods which have used in this paper to calculate the superfluid density, work 
in principle for many different systems. For example, one could use the methods used 
here to calculate the superfluid density in the system where another component of the 
two component Fermi gas occupies the first excited band [16], or in the case where an 
optical lattice is in a trap. These methods can also been used calculating the superfluid 
density of modulating gap states. 
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Figure 5. [Colour online] The BP/Sarma phase superfluid fraction. In figure (a) we 
show the superfluid fraction as a function of the polarization P = (n-f — n^)/(n-f — njj. 
Figure (b) shows the superfluid fraction divided by |A| 2 . Figure (c) shows the 
superfluid fraction as a function of the temperature. In figure (d) we show the 
superfluid fraction divided by | A | 2 . In figures (a) and (b) fc^T/J = 0.75 and in 
figures (c) and (d) P = 0.10. In all the figures all the hopping strengths are same i.e. 
Ja,u = Ja'.a' = J, the average filling fraction n av = 0.50, U/J= —6.0, and q = 0. 
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the FFLO superfluid fraction as a function of the total filling fraction, with three 
different interactions at zero temperature. In figure (b) we show p xx {n^ + njj as a 
function of the total filling fraction at zero temperature, when U/J = —6.0. In the 
both figures all the hopping strengths are same i.e. J a ^ a ~ Ja'.a' = J, q = IxX, and 
P w 0.30. The horizontal lines in the figures show the value 0. 
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Figure 7. Figure (a) shows the Sarma- and FFLO phase superfluid fraction as 
a function of the polarization at constant temperature. In figure (b) we show the 
Sarma- and FFLO phase superfluid fraction divided by |A| 2 . The parameters, which 
were used, are ksT/J w 0.23, U/J w —5.1, and the average filling fraction is 
n av = {n-f + n{)/2 = 0.2, All the hopping strengths were equal. Below P w 0.15 
q = and the phase is the Sarma phase and above P ss 0.15 q ^ and the phase is 
the one mode FFLO phase. The sudden drop in figure (b) at P w 0.5 in indicates that 
system becomes normal for higher polarizations. 
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Figure 8. [Color online] Phase diagrams with four different interaction strengths, 
when the average filling fraction n av — 0.2. The interaction strengths are from up left 
to bottom right:L//J « -3.7,17/ J « -4.4, U/J « -5.1, and U/J w -6.3. All the 
hopping strengths are equal. The color are as follows:BCS/Sarma=blue(black), stable 
FFLO=yellow(light grey), unstable FFLO=red(dark grey), normal gas=white. 
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approach to superfluidity of atoms in an optical lattice. J. Phys. B, 36:825, 2003. 
[27] R. Roth and K. Burnett. Superfluidity and interference pattern of ultracold bosons in optical 

lattices. Phys. Rev. A, 67:031602, 2003. 
[28] B. Wu and Q. Niu. Landau and dynamical instabilities of the superflow of Bose-Einstein 

condensates in optical lattices. Phys. Rev. A, 64:061603, 2001. 
[29] A. A. Burkov and A. Paramekanti. Stability of superflow for ultracold Fermions in optical lattices. 

Phys. Rev. Lett., 100:255301, 2008. 
[30] E. Taylor, A. Griffin, N. Fukushima, and Y. Ohashi. Pairing fluctuations and the superfluid 

density through the BCS-BEC crossover. Phys. Rev. A, 74:063626, 2006. 
[31] N. Fukushima, Y. Ohashi, E. Taylor, and A. Griffin. Superfluid density and condensate fraction 

in the BCS-BEC crossover regime at finite temperatures. Phys. Rev. A, 75:033609, 2007. 
[32] E. Taylor. The josephson relation for the superfluid density in the BCS-BEC crossover. Phys. 

Rev. B, 77:144521, 2008. 

[33] V. M. Stojanovic, W. V. Liu, and Y. B. Kim. Unconventional interaction between vortices in a 

polarized Fermi gas. Annals of Physics, 323:989, 2008. 
[34] M. Iskin and C. A. R. S de Melo. Superfluidity of p-wave and s-wave atomic Fermi gases in optical 

lattices. Phys. Rev. B, 72:224513, 2005. 
[35] I. Bloch, J. Dalibard, and W. Zwerger. Many-body physics with ultracold gases. Rev. Mod. Phys, 

80:885, 2008. 

[36] J. -P. Martikainen and H. T. C. Stoof. Spontaneous squeezing of a vortex in an optical lattice. 

Phys. Rev. A, 70:013604, 2004. 
[37] G. D. Mahan. Many-Particle Physics, 3rd ed. Kluwer Academic/Plenum Publishers, New York, 

2000. 



Superfluid- density of the ultra-cold Fermi gas in optical lattices 



21 



[38] E. M. Lifshitz and L. P. Pitaevskii. Statistical Physics, Part 2. Butterworth-Hcinomann, Oxford, 
2002. 

[39] A. L. Fetter and J. D. Walecka. Quantum theory of many-particle systems. Dover, New York, 
2003. 

[40] A. M. Clogston. Upper limit for the critical field in hard superconductors. Phys. Rev. Lett., 9:266 
- 267, f 962. 

[41] Meera M. Parish, Stefan K. Baur, Erich J. Mueller, and David A. Huse. Quasi-one-dimensional 

polarized Fermi supcrfiuids. Phys. Rev. Lett., 99:250403, 2007. 
[42] Y. Chen, Z. D. Wang, F. C. Zhang, and C. S. Ting. Exploring exotic superfluidity of polarized 

ultracold Fermions in optical lattices. 2007. 
[43] M. Iskin and C. J. Williams. Population imbalanced Fermions in harmonically trapped optical 

lattices. Phys. Rev. A, 78:011603(R), 2008. 
[44] D.-H. Kim, P. Torma, and J. -P. Martikaincn. Induced interactions for ultra-cold Fermi gases in 

optical lattices. Phys. Rev. Lett., 102:245301, 2009. 
[45] T. K. Koponen, T. Paananen, J.-P. Martikaincn, M. R. Bakhtiari, and P. Torma. FFLO state 

in 1, 2, and 3 dimensional optical lattices combined with a non- uniform background potential. 

New Journal of Physics, 10:045104, 2008. 
[46] J.-P. Martikainen, E. Lundh, and T. Paananen. Interband physics in an ultra-cold Fermi gas in 

an optical lattice. Phys. Rev. A, 78:023607, 2008. 



